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FOREWORD 


The research in this document is intended to contribute 
toward a quantitative assessment of the aerodynamic penalties 
on an aircraft with frost coated wings during takeoff. The 
frost problem is serious for both general aviation and air 
carriers. For air carriers it is an economic hardship because 
of the expense involved in removing frost prior to takeoff. 

For general aviation it is a potential safety hazard since 
takeoffs are permitted with frost. The objective of the UDRI 
research effort, of which this report constitutes the first 
task, is to quantify the safety hazards of frost so that 
realistic takeoff procedures and guidelines can be established. 

In the future, the model will be used for a twelve hour advance 
prediction of the density and thickness of morning accumulation 
of frost on an airfoil. The model will also be used as part 
of an aircraft simulation program under a frosted wing condition 
to determine the dissipation of frost during takeoff and climbout. 

This research was conducted by the UDRI for NASA/George 
C. Marshall Space Flight Center, Huntsville, Alabama, under 
the technical direction of Mr. Dennis Camp of the Space Science 
Laboratory. The support for this research was provided by 
Mr. John Enders of the Aeronautical Operating Systems Division, 
Office of Advanced Research and Technology, NASA/Headquarters. 
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SECTION 1 
INTRODUCTION 


In colder climates the overnight frost accumulation on an 
aircraft is a common occurrence. To ensure safe takeoff. 

Federal Air Regulations require that frost be removed from 
commercial aircraft. Frost removal is a costly and time con- 
suming process that perhaps could be avoided if better information 
were available about the frost formation, its dissipation, and its 
aerodynamic effects on takeoff. For general aviation, takeoffs 
are permitted with a frosted wing surface. Safety records 
indicate, however, that many takeoff accidents are attributable 
to the aerodynamic penalties associated with a frost roughened 
airfoil . 

Very little quantitative data exists concerning the aero- 
dynamic penalties associated with frost on an aircraft wing. In 
many cases a completely safe takeoff may be possible if the frost 
layer is sufficiently thin. In other cases where aerodynamic 
penalties are significant, a safe takeoff may still be possible 
at a reduced takeoff weight or with a sufficiently long runway. 

An additional area of concern for the aviation industry is the 
frost forecast. An accurate forecast of the frost severity on 
the eve of the frost would allow for adequate personnel, equip- 
ment, and supplies to be on hand for the frost removal excercise 
in the morning, or for some frost prevention procedure to be 
applied . 

The aerodynamic penalties resulting from frost accumulation 
have not been studied in detail due, apparently, to the three- 
dimensional nature of the flow problem. Also, the complexity of 
the frost formation process appears to be one of the reasons that 
a satisfactory numerical frost model is not in existence. Frost 
formation has been investigated experimentally for numerous flows 
over flat surfaces and deposition configurations by several 
investigators. However, little correlation exists between theory 
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and experiments and this is probably the reason that no general 
theory can be found. 

In this work, analysis of frost formation and its effect on 
takeoff aerodynamics proceeds by considering a more simple problem 
and then generalizing to more complex situations. The wing is 
first approximated as a flat plate and the flow around it is taken 
to be two-dimensional. The present study is made for three 
separate phases viz: 1) overnight formation and accumulation 

of frost on a flat plate; 2) modification of Phase I for the 
wing of the aircraft; and 3) the effect of accumulated frost 
on the wings during takeoff and the resulting aerodynamic 
penalties. Only the first phase is considered in this report. 

The study regarding the other two phases will be contained in 
future reports. The flat plate model developed in Phase I will 
be used for an airfoil in Phase II. It is anticipated that the 
two-dimensional approach can be extended to three dimensions once 
the complexities and the difficulties inherent in the present 
problem are analyzed. 

Frost is formed when moist air comes into contact with cold 
surfaces having temperatures below both the dew point and the 
freezing point. The determination of the overnight accumulation 
of frost on a flat plate is a complex heat and mass transfer 
problem. Radiation cooling to space off the top side of the 
metallic surface causes it to cool faster than the surrounding 
air. When the metallic surface reaches a frost point temperature, 
frost formation begins. Additional cooling then increases the 
frost density and its thickness. The frost accumulation eventually 
forms an insulation between the metallic surface and space and 
thus an equilibrium condition may be reached wherein the frost 
surface temperature remains constant. The important heat transfer 
processes shown in Figure 1 are the radiative cooling to space, 
conduction from the wing surface through the frost to the air, 
the natural convection of heat to the frost surface, and the 
internal enthalpy rate within the frost. The change of state 
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Air Velocity Profile 


x 



Figure 1. Schematic diagram of the heat transfer processes in 
the frost layer 



from water vapor to frost introduces an important heat source 
within the frost layer, the latent heat of sublimation. 

The nature of the frost formation process is sufficiently 
complicated so that it is difficult to predict the rate of 
formation and its density at a particular time. Many of the 
important characteristics of the frost depend on how it was 
formed. This history dependent nature of frost formation forbids 
the use of common types of correlations for predicting heat 
transfer rates used in other heat and mass transfer problems. 
Ordinarily, if the surface temperature is known, the heat trans- 
fer from the moist air to the surface can be calculated. Thus, 
if it is possible to describe adequately the heat transfer 
through the frost, then one can find the frost surface temperature 
by matching the heat transfer rate through the frost with the 
convective heat transfer. The description of the heat transfer 
rate through the frost requires a knowledge of the thermal con- 
ductivity and the history of formation of the frost which have 
not been adequately treated to date. Section 2 of this report 
is devoted to a treatment of the thermal conductivity and a new 
expression is presented. Section 3 analyzes the heat and mass 
transfer coefficients and gives modified relations which can be 
used more accurately for the transfer coefficients at the frost 
surface. In Section 4, a model is presented for predicting the 
growth rate of a frost layer using the analysis of the previous 
sections. Since the equations so developed cannot be solved 
analytically, numerical techniques have been applied to solve 
them. These techniques are also presented in Section 4. In 
Section 5 the results from the numerical model are compared with 
the experimental data available on frost formation. Considering 
the complexity of the problem, the results are very encouraging. 
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SECTION 2 

DERIVATION OF THE FROST THERMAL CONDUCTIVITY 


The thermal conductivity of the frost layer plays an important 
part in its structure and rate of formation. A number of authors 
addressed the problem of computing the frost thermal conductivity 
(References 1, 2, 3, 4, and 5). In this section, the various 
approaches used by these authors are examined by analyzing the 
underlying assumptions of each treatment. To fully understand 
these assumptions it is necessary to begin with a discussion of 
all possible heat transfer processes 'within the frost layer and 
to determine which processes are significant and which can be 
safely neglected. From this perspective, the different approaches 
taken by the authors can be evaluated and their results compared 
with experimental data. Furthermore, the range of environmental 
conditions will be determined for which a particular approach is 
realistic and the limitations of each approach can then be 
deduced. It will be shown that none of the present approaches 
are sufficient for a general frost model. As a result, a new, 
more comprehensive method of calculation the frost thermal 
conductivity based on the experimental data will be developed. 


2.1 THE HEAT TRANSFER PROCESSES WITHIN THE FROST LAYER 


The various approaches to modeling of the frost temperature 
distribution and thus also the frost thermal conductivity can 
be derived from expressing the different heat transfer processes 
within the frost layer as shown in Figure 1. The expression 
for the energy equation for a control volume is given by 
(Reference 6) 


pf §§ = §i + ft + 7 - (k7T> + - 7 ^ 


( 1 ) 


where 


is the total derivative =* + U ' 

Dt 3t a 9x 

p f is the mass density of the frost layer (g/cc) , 

t is time (s) , 
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i is enthalpy (per unit mass) , 

P is pressure (N/m 2 ), 

Q is internal heat generation within the layer (w/m 3 ), 

k is the thermal conductivity of ice/air 

structure in the layer (w/m°C) , 

T is temperature (°K), and 

q r is the radiation heat flux vector (w/m 2 s) . 

For heat transfer within a frost layer the following 
assumptions are made. 


1 Within the frost layer, the temperature and 
the pressure are in a quasi-steady state, 
i.e., within a time interval At 


2 


p f C pf 


3T 

3t 


DP 

Dt 


0 


( 2 ) 


This assumption is based on experimental 
data by J. White (Reference 7) . White has 
shown that the temperature and the pressure 
in the frost layer are at most slowly 
varying functions of time, partly due to 
the isothermal conditions of the wall and 
the upper limit of the frost surface temp- 
erature at melting point. Then, the energy 
storage rate is small compared with the 
heat flux; thus the justification of the 
assumption of quasi-steady state. 

The internal heat generation rate produced 
by the phase change is given by 


30 _ T 3 *d 

3t L 3x 


where L 


L 


e 



= L or L and 
e s 

is the latent heat of water 
evaporation (J/g) , 

is the latent heat of ice 
sublimation (J/g) , 

is the water vapor mass flux 

within the frost layer (g/m 2 s) , 


(3) 


and 


x is the distance from the wall (m) . 


3 


The enthalpy change produced by forced air flow 
through the frost structure is given by 


GC 
a p 


3T 

9x 


p .U 
d a 


3i 

3x 


( 4 ) 
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4 


5 


where G is the air mass flow rate per unit 
a area (g/sm 2 ) 

C is the specific heat of air ( J/g°C ) 
P 

U is the air velocity (m/s) . 
a 

The heat conduction is one-dimensional 
(through the layer to the wall) and the 
effective thermal conductivity, k , is a 
function of the combined heat conducting 
proportions of air and ice. This gives 

< kVT > = (k e -Hr >• 

The radiation heat flux is one-dimensdonal , 



where k = radiation thermal conductivity 

1 based on the Stef an-Boltzman law 
and the geometric view factors. 


and 


(5) 

( 6 ) 


Under these assumptions the energy Equation (1) becomes 

driu 


dx ^ ^ k e + ^r^ dx^ J " dx '"'a w p dx 
A familiarity with the above energy equation by order of 
magnitude calculations is needed. More specifically, four 
terms, k , k , m,. and G , will be investigated in the above 
equation. It will be demonstrated that the heat flux represented 
by the radiation and the forced air enthalpy transport are 
negligible in comparison to the thermal conductivity of air- 
ice structure and the latent heat release of the waper vapor. 


dT 


■] = -L 


- G C 


dT 


(7) 


Since it will be shown that the thermal conductivity of air- 
ice structure contributes the largest of the four heat fluxes, 
the extreme values of k Q are used as a reference comparison for 
the order of magnitude calculations. The extreme values are 
the thermal conductivity of ice (Reference 8) 


k i = 630/T 


and the thermal conductivity of air (Reference 6) 

m 1/2 


k = 2.646 X 10 
a 


-3 


1 + ( i45) 10 - 12 /T 


(8) 

( 9 ) 
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2.1.1 Air - Ice Thermal Co nduc t i y i ty , k ^ 

As the frost contains air and the crystals of ice, 
the conductivity of frost should be somewhere between the thermal 
conductivities of air and ice. In the limit, depending on its 
density, it should approach the thermal conductivity of air or 
ice. Thus, the thermal conductivity of frost is defined in terms 
of the weighted functional relationhip between the thermal 
conductivities of air and ice based on the density and the struc- 
ture of the frost. Biguria and Wenzel (Reference 3 ) have compiled 
several theoretical models for different systems to formulate the 
effective thermal conductivity of the frost based on various 
assumptions for the frost structure. These expressions as applied 
to the air and ice phases are given as follows. 


1 . 


Resistance in series for minimum possible 
conductivity : 


1 

'em in 



/ 


where 3 is the porosity of frost. 


2 . 


Resistances in parallel for maximum possible 
conductivity : 


emax 


(1-3') k. + 3 k - 

_L d 


( 10 ) 


(ID 


3 . Russel equation for porous media where the 
solid (ice) is in a continuous structure and 
there is a distribution of cubical pores 
arranged in a simple cubic lattice: 


4 . 


^e ^ a3 2/3 + 1 - 3 2/3 

k i a (3 2/3 - 3 ) + 1 -3 2/3 +3 

with a = k /k. and 0 = 1-3 ■ 
a' x 


Maxwell - Rayleigh equation for the case 
of fluid pores (air) distributed in 
a continuous solid (ice) : 


( 12 ) 



( tt 




(13) 
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5 . 


Maxwell- Rayleigh equation for the 
case of solid pores (ice) distributed 
in a continuous fluid (air) : 



6 . 


If one phase of the constituents (say 
ice or air) is not spatially continuous 
the Brailsford-Ma jor equation gives: 



+ 


+ 


| (36 - l)k a + 


x 


(3 6- 1) k + 

ci 


8k k 

a 




(30 - 1) k ± 

(30 - 1) k ± ) 2 


7. Woodside equation for a cubic lattice of 
uniform solid spherical particles (ice) 
in a gas (air) : 



i 4 ) 1/2 

where a = 

for 0< G < 0.5236. 


( 14 ) 


(15) 


(16) 


In their observations of frost formation. Brain, et al . 

(Reference 1) found that the initial frost dendrites are: spherical 
in shape at about 5 to lOy in diameter. As smooth frost forms, 
the diameters become about 20 to 50 y and the ice dendrites begin 
to mesh together. Biguria and Wenzel (Reference 3 ) observed 
that initial frost was rough, consisting of ice trees and air 
spaces. They assumed that parallel heat transfer could be 
dominant up to a frost density of 0.02 g/cc. Then from 0.02 g/cc 
to 0.05 g/cc, the thermal conductivity was observed to decrease 
since parallel heat transfer was no longer valid when the frost 
formed a close-knit mesh of dendrites. Then at densities 
greater than 0.05 g/cc, the dendrites begin to enclose air pockets. 
Thus a realistic frost structure model should be spatially 
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continuous both in the air and the ice phases. Note that none 
of these seven theoretical thermal conductivity expressions 
directly provide for such a structure. 

Given the complexity of deriving an expression of 

k for a realistic frost structure, one wonders if it is useless 
e 

to list Equations (10) to (16) and attempt a theoretical approach 
toward the frost thermal conductivity. Later, an entirely 
empirical approach will be found to be of little value also. 
Finally, a semi-theoretical approach which makes indirect use of 
Equations (10), (11), (13), (14), and (15) will be derived to reflect 
a realistic frost structure. At this point it is important to 
understand that k Q is the air-ice thermal conductivity. The 
frost thermal conductivity is an expression that will be derived 
to take into account the other heat flux terms in Equation (7) . 

2.1.2 Radiation Effective Conductivity, k r 

The radiation effective conductivity can be shown to 
be negligible for the size of the ice crystals and the temperatures 
in the frost layer by the following argument. The radiation 
effective conductivity as given by Laubitz (Reference 9) is 

k r = 4 ctT 3 £ | (1 - 0 2/3 + 0 4/3 ) , (17) 

— 8 2 4 

where a is Stefan-Boltzmann constant - 0.56697 X 10 w/m /°C , 

£ is emissivity = 0.985, 

b is the linear dimension of ice crystals, and 
0 is the fractional volume of ice fragments - 1-3. 

Substituting in the above expression, the typical 
values p f = 0.13 g/cc and T - 266 °K (Reference 10) we obtain 

k < 1.43 x 10“ 4 w/m°C . 
r — 

However, from Equations (8) and (9) at the same temperature we 
have , 

k^ - 2.37 w/m°C and k a = 0.0236 w/m°C. 
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This indicates that the radiation effective conduc- 
tivity is some two orders of magnitude less than the thermal 
conductivity of air. In addition, a typical experimental data 
of frost thermal conductivity appears to show a noise level 

— 3 

around 10 w/in 0 c or more. Therefore the radiation effective 
conductivity is considered negligible. 

2.1.3 Water Vapor Effective Conductivity, k 

The concept of the water vapor effective conductivity 

drti-, 

is obtained by assuming the energy term, L , in Equation (7) 

obeys the diffusion equation and meets the condition of water 
vapor saturation in the frost layer. This concept is, of course, 
invalid when the frost layer becomes supersaturated or sub- 
saturated. Whether these conditions exist or not depends on the 
significance of nonequilibrium dynamics versus a strong tendency 
toward an equilibrium state within the frost layer. Since in 
frost formation whereby the water vapor flux enters the frost 
layer, the state of subsaturation is quite unlikely. If the 
water vapor flux into the frost layer is so rapid that homo- 
geneous nucleation occurs, then a state of supersaturation 
exists (Reference 11) . But because there are several nucleation 
sites within the layer to prevent homogeneous nucleation, one 
expects there exists a critical wall temperature above which 
a state of supersaturation is quite unlikely. Experimental 
evidences of these observations will be shown later. 


If the frost layer can be assumed to be in the 
saturated state, then the water vapor mass flux is given by the 
following diffusion equation for the frost (Reference 5 ) . 



D 3 

(1-X)t s 



(18) 


where 


D = 1.198 X 10~ 5 Tl,75 ( p atm / p ) is the diffusion 
coefficient (Reference 12 data) , 


X is the relative concentration (moles f^O/moles air) , 
3 is the porosity = ( “ Pf)/( p i “ p a ^ # 
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t is the tortuosity, and 

P v is the saturated water vapor density (g/cc) . 

The porosity accounts for the decreased effective 
cross sectional area for diffusion and the tortuosity, generally 
taken as 1.1 for frost, accounts for the increased path length 
the molecules must travel. The assumption (whose verification is 
apparent later in Figures 2 and 3) that the water vapor at the 
frost surface is saturated implies that the water vapor mass flux 
can be made to follow the temperature gradient through the gas law. 


P , = p R T , 
v v v 

and the Clapyron Equation, 


P = P 
v v 


exp 




R T 
v 


R V T 


] 


(19) 


( 20 ) 


where 


p v 

T* 


V 


Setting L = L e , P 
P 

v 


is the latent heat of ice sublimation or 
water evaporation (J/g) , 

is the referenced pressure (N/m 2 ), 

is the referenced temperature (°K), 

is the water vapor gas constant = 4.6150 X 10 6 erg/C 2 g, 
is the water vapor density (g/cc) , and 

is the water vapor pressure (N/m 2 ). 

* 

s ' P v 

= 610.7 exp (22.4959 - 


= 610.7 N/m 2 , and T* = 273 °K, we obtain, 

L_ 


R T 
v 


-) 


By differentiating the gas law equation, we obtain, 


3P, 


dP 


v 


dT 


( 21 ) 


3x 


R T 
v 


dT 


R T‘ 
v 


dx 


The Clapyron equation for P^ is differentiated with 
respect to temperature and is substituted along with Equation (21) 
into Equation (18) to get, 


m 


D6 


TT-xTtT 



dT 


dx 


( 22 ) 
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Figure 3. Water mass flux versus distance in frost layer 
(Yamakawa, et al. data. Reference 13) 
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Since the water vapor mass flux is now directly 
related to the temperature gradient, a thermal conductivity due 
to the water vapor latent heat flux is defined by 

nuL L D$ / P \ /L \ 

k v 5 (dT/dx) = (1-X)T S " J <23) 

Equation (23) is seen to be closely linked to the 

saturation conditions through the term (dP,VdT)/R T in Equation (21) 

2 ^ 

The second term in Equation (21) , p^/R^T , was found to be quite 
small in comparison. If supersaturation exists, then the 
Clapyron equation is no longer valid' and a new equation for P v 
would have to be derived for this state. Fortunately, it was 
found this is not necessary, as the following experimental 
evidence will show. 


According to Equation (22) , if experimental values 
of the frost density and temperature distribution are known, 
then the water vapor mass flux can be calculated as a function 
of the "x" variable as shown in Figure 1 within the frost layer. 
Observations by several authors (References 1, 2, 3, 13, and 14) 
indicate that the frost density is nearly spatially invariant in 
the "x" direction. This means the total water mass flux, 
consisting of water vapor and nucleated water drops, is given 
approximately by, 


dm, 
t 

dx 



x 


m 


ex£ 

x 


s 


(24) 


where i^ eX p is the experimentally calculated value of the water 

mass flux into the frost surface from the surrounding air and 

x is the frost thickness. If one observes that m. < m. for all 
s t — d 

values of x, then no nucleated drops have formed? this means 

a supersaturated state is unlikely. If one observes that 

• • 

m t > m d ^ or some values of x, then nucleation has occurred, and 
thus supersaturation might be possible. Note that homogeneous 
nucleation and nucleation on nucleating sites in the frost layer 
cannot be experimentally distinguished. Thus we cannot state 
definitely if supersaturation has occurred. If m^ = m^ at the 
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frost surface or at distance, x , then we have a good method for 

s 

predicting the frost density growth rate for which a formula 
will be derived in Section 4 . 

Figure 2 shows a plot of m^ and m t calculated from 

an experimental run of Brian, et al . (Reference 1), where the 

wall temperature was at 80 °K. The steep rise of the m^ curve 

is due to the temperature dependent relationship of D and P 

in Equation (22) for the temperature increase from 80 °K to about 

265°K. In comparison to the m t curve it is probable that some 

supersaturation has occurred, given the magnitude of the difference 

between m. and flu . At distance x from the wall we note that 
t d s 

m^ is equal to m t . A data set with a more reasonable temperature 
range applicable for frost formation on an aircraft can be obtained 
from Yamakawa, et al . (Reference 13). Here the range of tempera- 
ture, in one specific case, is from -22°C to -3.3°C. Although 
the experimental temperature distribution within the frost is not 
available, indications are that for this small temperature range, 
the temperature profile can be roughly approximated as a linear 
function of x. Thus, the temperature gradient for Equation (22) 
is given by 

~ = 66 . 78 o K/cm (25) 

dx x 

s 

as obtained from experimental data in Reference 13; where x g = 

0.28 cm. At this frost thickness the frost density is 0.1110 
g/cc and the ambient absolute humidity is 0.0049 as obtained from 
the data. Substituting these values into Equation (22) gives the 
m^ curve shown in Figure 3. For the top half of the frost layer/ 
the curve agrees closely with the m t curve calculated from 
experimental data. While for the lower half of the frost layer, 
m^ is greater than m^. These observations mean that down to a wall 
temperature of -22 °C we can confidently say the frost layer is in 
a saturated state. In addition, at the distance x g , we find that 
m^ - m^. The conclusion is that the water vapor thermal conduc- 
tivity is valid for most frost formation situations and an accurate 
method for calculating the water mass flux entering the frost 
surface has been obtained. 
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Now k^ can be compared directly with k^ and k^ for 
an order of magnitude analysis. From Reference 10, with the same 
data used in evaluating the radiation effective conductivity, 
we obtain at the frost surface, with equation (23)., 

k * 0.0111 w/m°C • (26) 

v 

But since k = 0.0236 w/m°C, the water vapor thermal conductivity 
a 

cannot be ignored at low frost density. At higher frost density, 
k v actually decreases due to the porosity term in Equation (23) . 

At close to ice density the term k^ = 2.37 w/m e, C indicates the 
dominating influence of k Q . So far, the order of magnitude 
calculations show that particular attention must be devoted to 
developing the air-ice thermal conductivity which would include 
a complicated frost structure modeling and perhaps also the water 
vapor thermal conductivity at low frost density. 

3T 

2.1.4 Forced-Air Enthalpy Rate Term, G^ C ^ 

The air mass flux, G , is difficult to determine 

a, 

because it is strongly dependent on the frost structure. If 
the frost structure consists primarily of ice cylinders that 
penetrate deeply into the ambient flow of air, then perhaps G a 
can be approximated conservatively to the upper limit, e.g. , 
p V*, . But that situation is unlikely since, as shown above, 
the frost structure is not simple. That is, the initial frost 
may form ice trees, but the frost thickness is so thin that it 
barely penetrates the momentum boundary layer. This means the 
velocity, v , is some small fraction of the free stream velocity. 

As the frost thickens, the frost density will also increase, 
causing a close knit mesh of ice dendrites. This frost structure 
would eliminate any forced-air flow through the frost. It is 
noted the wall is impermeable, so there is no suction or blowing 
underneath the frost layer. In the study of heat and mass 
transfer coefficients in Section 3, for the forced convection 
case, we show that the effective frost surface area for the 
turbulent heat transfer coefficient remains a constant for all 
values of the frost density and thickness. But for the laminar. 
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natural convection case, there is no frost surface effect on 
the heat transfer coefficient. Thus it is concluded that the 
momentum boundary layer begins at the frost surface for the 
laminar natural convection case and begins at some constant 
minute distance below the frost surface in the turbulent 
forced convection case. In this situation, the only mass flux 
within the frost layer is the total water flux, m^, which is 
set equal to G . 

cl 

Consequently from Brian data, as also used to evaluate 
k and k r , a high estimate is made of the forced-air enthalpy rate 
term in Equation (7) for comparison with the fh^ term by using 
the measured temperature gradient for the water vapor mass flux 
at the frost surface. The result is, 

= 8.44 X 10“ 4 w/cc . (27) 

s 

This can be compared to a low estimate of the latent 
heat contribution calculated by setting 



= 2.91 X 10 2 w/cc • (29) 

Comparable results were also obtained within the 
frost layer. Therefore the heat transfer rate by the forced- 
air enthalpy rate is much lower than the latent heat release 
rate within the frost layer. In fact it also appears that the 
heat transfer by radiation within the frost layer is of the same 
order of magnitude as the heat transfer by an enthalpy change 
due to the water vapor flux. Since the experimental data is not 
accurate to four significant digits as would be required both 
by the forced-air enthalpy rate and the radiation heat rate, the 


dx 


in Equation (7) to get 


m 


d L 
s s 


. „ dT 

d a dx 
s 
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terms for the conductive heat rate and the latent heat release 
rate are only retained in Equation (7) The result is 


( 30 ) 


[ k arl - . L dlh a . 

dx [ K e dx J ^ dx 

An expression for a frost thermal conductivity can 
be derived if Equation (2.3) is used and substituted in the above 
equation and integrated. Thus, 


with K * k + k v , - (32) 

where q Q is a constant heat flux at the wall and K is the thermal 
conductivity of frost. With the simplified equations above, 
all other approaches obtained from the literature for calculating 
the frost thermal conductivity can be explained. 


2 . 2 APPROACHES FOR CALCULATING K 
2.2.1 Brian, et al . Approach 

If the heat flux and the temperature gradient are 
measured, as was done in the frost experiments by Shah (Reference 2) 
and Brian, et al . (Reference 1), a frost thermal conductivity 
can be calculated easily from the above equations. It is 
important that the heat flux should be measured at the wall as 
is done by Shah (Reference 2) as well as by Brian, et al . 

(Reference l) rather than at the frost surface. The reason is 
because, due to the latent heat contribution, the heat flux 
becomes lower at the frost surface than at the wall. 

Some empirical expressions for the thermal 
conductivity of frost, based on the experimental data by Shah 
(Reference 2), Brazinsky (Reference 15), and Brian, et al . 
(Reference 1), as developed by Brian, et al. ( Reference 10) in 
SI units are 
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2.405 


(33) 


K 1 = 


where P 


K 2 = 3.595 


, a -5 m 1.272 . no , , n -5 m 1.74 
x 10 T + 3.931 x 10 p f T 

< 0.025 g/cc and T < 255° K; 

-15 5 44 -12 4 84 

x 10 + 2.042 x 10 (P f -0.025)T * , 


(34) 


where p^ > 0.025 g/cc and T _> 255° K; and 

-5 1 441 -8 3 055 

K 3 = 1.564 x 10 °T * + 4.252 x 10 ° (P f -0 . 025) T J * UDO , (35) 


where P^ > 0.025 g/cc and T < 255° K and 

where P^ stands for the average density of the frost. It may be 
noted that this empirically developed thermal conductivity is a 
linear function of the frost density. Later on it will be shown 
that this relationship is valid only within the experimental 
range of frost density up to 0.13 g/cc. 

Serious problems are encountered if one attempts to 
use the above type expression for thermal conductivity with other 
data. This is because a number of factors such as the thermal 
conductivities of water vapor, air, ice, and the frost structures, 
which play important parts in the expression for the thermal 
conductivity of the frost, have not been considered in the above 
model . 


2.2.2 White * s Approach 

Another approach to determine the frost thermal conduc- 
tivity developed by White (Reference 4 ) using the empirical 
relations to fit the Shah data (Reference 2 ) where k v of Equation 
32 is given by 


k 

e 


k 


v 


D 


ef f 


R 



and 


(0.1684 + 35.39 P f /P water 


X (-0.01356 + 0.0001579T), 


(36) 

(37) 


and 
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2 

where D e ff = 100 111111 / sec - Although this, approach is more 
theoretical than that of Brian, et al. , since it takes into 
account more variables, it is still hampered by the fact that k v 
and k^ are treated as empirical terms. Another problem is that 
this approach is valid only in the experimental range of Shah's 
data (Reference 2) which is 0.13 g/cc. 

2.2.3 Biguria and Wenzel's Approach 

Since the water vapor diffusion occurs only in the air 
portion of the frost, Biguria and Wenzel (Reference 3 ) suggest 
that if the effective air thermal conductivity instead of the true 
air thermal conductivity is used, one can expect better results. 
The effective air thermal conductivity can be obtained from the 
relation 


E>P / 

k . = k + - — ( 
eff air a P - P I 


L P 

S V 


R V T 



k + 
a 


t k 
s v 


6 


(38) 


which is based on Equation (21) of Reference 3 and Equation (23) . 

It may be noted that for obtaining the air effective conductivity, 

the tortuosity and the porosity expressions in the vapor effective 

conductivity are neglected as they have no use as corrections in 

the air portion of the frost. According to Biguria and Wenzel 

(Reference 3) one can obtain a better expression for k by using 

k nr . instead of k in Equations (10) to (16) . For frost 
eii air a 

densities greater than 0.05 g/cc, Biguria and Wenzel claimed that 
a good fit to their experimental data was obtained by using this 
approach and Equations (15) or (16) . The experimental thermal 
conductivity of frost was obtained, however, by measuring the heat 
flux at the wall, the thickness of the frost, and the wall and 
frost surface temperatures. In order to apply the theoretical 
thermal conductivity equations to the data, Biguria and Wenzel 
implicitly assumed the frost layer has uniform temperature. 
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structural , and density distributions. These assumptions require 
the wall temperature to be a high temperature at around 250° K, 
the frost to be formed in a specified structure, and the ambient 
absolute humidity to be in a specified range near saturation. 

In contrast, the data of Brian, et al. (Reference 1) and 
Shah (Reference 2) typically has the wall temperature at 80° K, 
more spherical ice formations than ice trees , and an ambient 
absolute humidity at a fraction of the saturation level. Thus it 
is expected the approach by Biguria and Wenzel will not fit 
satisfactorily to the data of Brian, et al. (Reference 1 ) and 
Shah (Reference 2 ) . This actually turned out to be the case. 

The basic disadvantage in this approach is the requirement of a. 
uniform temperature and structural distribution of the frost. 


2.2.4 Jones and Parker's Approach 


By assuming that the frost density is spatially 
invariant so that the amount of water vapor being frozen is the 
same at all locations in the frost layer, the change in the vapor 
mass flux with distance is 

d *d _ (28) 

"dx x~ 
s 


where m d s is the water vapor flux driven by the temperature 
gradient at the frost surface. Jones and Parker (Reference 5) 
have used this relationship to model the heat transfer mechanism 
in the frost layer, and thus also the model of the frost thermal 
conductivity. Substituting the above expression into Equation (30) 
and integrating, we get. 


k 


e 


dT 

dx 



x + q 
^o 


(39) 


where q Q = h H (T a - T g ) + h m L s (oi a - u> g ) (w/m 2 ), 

h H is the heat transfer coefficient (w/m 2 °C) , 

h is the mass transfer coefficient (g/m 2 sec°C) , 
m 

T is the ambient air temperature (°K), 

u 

T g is the frost surface temperature (°K), 
oj is the absolute ambient humidity, and 
03 g is the frost surface humidity. 
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It may be noted here, k^ , which is computed in the 
usual way as in Equations (10) to (16) , does not contain the 
expression for k^^ Also, Jones and Parker (Reference 5) 


air 


used the incorrect expression for k g by using Brian, et al . ' s 
(Reference 1) empirical frost thermal conductivity. This 
approach has two difficulties. If a linear vapor mass flux with 
distance which releases the latent heat is assumed, then, 
because of what has been shown in Figures 2 and 3, a supersaturated 
state would have to exist in Brian, et al.'s data (Reference 1) 
and a subsaturated state would have to exist in Yamakawa's data 
(Reference 13) within the frost layer. This dual state is 
physically unlikely, especially if several ice crystals also exist 
along side of the moist air within the frost layer. The second 
difficulty is that this approach cannot be compared with experi- 
mental data unless the term m, can be calculated from each data 

d 

s 

set. Usually the values of ifi^ are absent in the literature 
available. s 

2.2.5 Summary of the Four Approaches 


If one examines carefully the aforementioned approaches 
which have been established on the basis of a controlled environ- 
ment, one will find that they are not satisfactory for the air- 
craft frost formation problem where the environmental variables 
have different ranges. The validity of the first approach, given 
by Brian, et al. (Reference 1), is restricted to frost densities 
less than 0.13 g/cc and very low wall temperatures and ambient 
humidities, as indicated in Table 1. If the empirical frost thermal 
conductivities based on ice density are extrapolated, the sharp 
disagreement with the actual thermal conductivity of ice is 
apparent in Figure 4. Similarly, if the same is done on the 
basis of air density, the resulting curves in Figure 5 show 
the disagreement with the actual thermal conductivity of the air. 

It may be noted that the effective conductivity of water vapor, 
which becomes negligible for temperatures less than 225° K 
due to the rapid decrease in D and P v , was not considered in the 
thermal conductivity of air. On the basis of the same arguments. 
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TABLE I 


Summary of Approaches to Calculating 
Frost Thermal Conductivity 


Range of Modeling Modeling of the 


Approach 


Application 

Technique 

Frost Structure 

Brian, et al. 
(Reference 1) 

p f < .I3g/cc 

T » 80° K 
w 

Empirical 

None 

White 

(Reference 

4) 

p f £ . 13g/cc 

T = 80° K 
w 

Semi-empirical 

Vapor diffusion 
in frost layer is 
postulated 

Biguria & 
Wengel 


.05g/cc<p f <P ice 

Theoretical 

Simple Frost 
structures is 

(Reference 

3) 

T > 250 ° K 
w 


postulated 

Jones & 
Parker 
(Reference 

5) 

Pf£ ,13g/cc 

Semi-empirical 

Vapor diffusion in 
frost layer is 
postulated 

UDRI 

Approach 


p a < p f < p ice 

80° K<t <265° K 
w 

Mostly theoretical 
partly empirical 

Complicated frost 
structure is postulated 
for vapor diffusion, 
geometrical shapes of 
ice dendrites, and for 
frost aging. 


Limitations 


Average Thermal 
conductivity will not 
fit data at T other wall 
temperatures. 

Same as Brian , et al. 
(Reference 1) 


Same as Brian# et al. 
(Reference 1) 


Super saturation within 
frost layers is un- 
likely. 



Thermal Conductivity, (w/m °K) 





Thermal Conductivity, (w/m °K) 



Temperature °K 


Figure 5: Thermal conductivity versus temperature at air density 


the approach given by White (Reference 4) is found unsatisfactory 
for our problem. Since the state of super- or subsaturation 
within the frost layer is unlikely to exist, Jones and Parker’s 
approach (Reference 5) is not being used in our analysis. When 
the available experimental data could not be fitted to the 
theoretical approach suggested by Biguria and Wenzel (Reference 3) 
an attempt was made to see if the experimental values of frost 
thermal conductivity given by Brian, et al. (Reference 1) could 
lie between the curves represented respectively by Equations (13) 
and (14) for either spherical air pores or spherical ice particles 
The encouraging results in Figure 6 showed that the frost thermal 
conductivity is a linear function of porosity or frost density 
for porosities greater than 0.85. This gave a motivation to 
propose a new model, based on Biguria and Wenzel's theoretical 
approach (Reference 3) , but which includes frost structure para- 
meters which could be empirically derived to fit Brian, et al.'s 
data (Reference 1) , as well as other data. 

2.2.6 The UDRI Approach 

The proposed model makes the following assumptions 
about the frost structure as shown in Figure 7. At low frost 
density, or at high porosity, two types of frost structure pre- 
dominate. One is the ice cylinders created by the internal 
diffusion of water onto the ice, which result in a parallel con- 
ductive heat transfer. The other portion is the ice spheres 
created by nucleation of water vapor, resulting in a much lower 
conductive heat transfer. The total structure of the frost is 
then the random mixture of ice cylinders and ice spheres 
(Figure 7a). At high frost density, or low porosities, completely 
different dual structures begin to take shape. In contrast to the 
low density case, spherical air voids are formed in place of ice 
cylinders (Figure 7b) . This results in enhanced thermal conduc- 
tion. Also, in place of the ice spheres, stratified layers are 
formed. The total frost structure is then a random mixture of 
air bubbles and ice layers. 
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Total Thermal Conductivity (w/m °K) 



Figure 6 : Thermal conductivity versus frost 

porosity at 211 °K 
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Random mixture of ice planes and air bubbles 
at low porosities or high frost densities 

(b) 


Figure 7. Frost Stucture Model of the present work 
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With such a model of the frost structure, the equations 
for air-ice thermal conductivity presented earlier are combined 
in two ways. First, an attempt was made to arrive at the upper 
limit and the lower limit of thermal conductivity as a function 
of porosity. The second part attempts to ascertain the combined 
contribution of the upper limit and the lower limit of the thermal 
conductivity to the frost thermal conductivity as a function of 
temperature. Noting that the thermal conductivities for ice 
cylinders and for air bubbles are close to each other for all 
porosities, a simple interpolation rule is used to obtain the 
upper limit of the thermal conductivity. 


k u - (1-B) + 6k c (upper limit) , 

k b = k, [l - 26 (Hv)]/ l 1 + 6 (m)] 


(40) 

(air bubbles) , (13) 


= (1 - 3)k i + 3k ff air (ice cylinders) , 


( 11 ) 


a = k « . /k . 

eff air l 

Likewise, the lower limit of thermal conductivity is 
formed by an interpolation between thermal conductivities for 
ice spheres and ice planes. 

k, = (1 - B)k + 3k (lower limit) (41) 

1 p s 

k s = k i J3 + 2 3 ( ot - 1 )J / [ 3-3 ( a ~ -) j(ice spheres) , (14) 

k_ = ^i ^eff air (ice planes) . (10) 

p (1 - 3) k ~ : + k. 3 

eff air 1 

To determine the combined contribution of the upper 
and the lower limits of the thermal conductivity to the frost 
thermal conductivity, the random mixture model of Brails ford 
and Major given by Equation (15) is utilized. The present 
model of thermal conductivity of frost given as 


K = 


(33 - 1) k. + (3© - 1) k + \( 

c 1 c u |_V 


1/4 {( 

+ (30 C - 1} k u) 2 + 8k l k u] 1/2 } 


( 3 3 - 1) k. 

c 1 


(42) 
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is used to combine k, , k , and 0 . The proportion of the frost 

Ji. u c 

volume, B c , representing ice spheres and ice planes given as a 
parabolic function of the porosity B is 

3 = a + bB + cB^ , 


( 43 ) 


where a, b, and c can be constant or functions of temperature. 

The other portion of the frost volume representing ice cylinders 
and air bubbles is given by 

0 = 1-6 , < 44 > 
c c 

It is plausible to assume the frost takes on a completely 

spherical air voids structure when the porosity in the frost 

approaches zero or the density of the frost approaches that 

of ice. This assumption when translated into boundary condition 

for Equation (43) , gives 3 =0 for 3 = 0 and thereby 'a' is 

c 

found to be zero. Another plausible assumption is that all 
frost structures converge to the same structure as the porosity 
approaches zero regardless of the frost temperature. It gives 
'b' as a constant and ' c' will then be a strong function of 
temperature. Since each curve in Figure 8 is at a different 
temperature, it is possible to determine the values of ' b' and 
' c' from it. The result indicates that ’b’ is indeed a constant 
and 'c 1 is a function of local temperature of the frost, hence the 
value of 


3 becomes 
c 


B c = 6 


jo. 9 + 0 . 


1 + 


80 2 - 


(357) 


- 3 |- 


(45) 


As shown in Figure 8, it gives a multiple correlation 
of 0.99 and a standard deviation to the third decimal place 
as a fit to Brian, et al . ' s data (Reference 1) . In order to 
compare the frost conductivity, as predicted from the proposed 
model, to other experimental data, an average frost thermal con- 
ductivity is needed. The average frost thermal conductivity, 
denoted by K, is defined as 


f s 


K = 


'w 


K dT 


(46) 


T - 
s 


w 
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Total Thermal Conductivity, (w/m°K) 
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Figure 8. Comparison of the Present Frost Thermal 
Conductivity Model with Experimental 
Data of Brian, et al . (Reference 1) 



where T g is the frost surface temperature and T is the wall tem- 
perature. The calculated value of K agreed excellently with 
Brian, et al . ' s data (Reference 1). But Yamakawa, et al . 

(Reference 13) and Nakamura T s data (Reference 14) gave biased 
values on the low side for K. The reasons for the discrepancy 
seem to be the following. The effects of frost aging and the size 
of the water droplets arriving at the frost surface cause the 
formation of various types of frost structures. In Brian, et al . 
data (Reference 1) , where 

T t _ = 80° K and $ - 20%, (47) 

w a 

the extremely low temperature coupled with low relative humidities 
is found to cause the formation of several large size water drop- 
lets, which may freeze instantaneously on the wall. This con- 
clusion was derived from Rosner and Epstein (Reference 11) in the 
case of fog formation near cold surfaces. As the frost surface 
temperature increases with time, the smaller water droplets arrive 
at the frost surface and due to sharp temperature gradient, nucleation 
takes place within the frost layer. Also, as the frost density 
increases, the ice dendrites begin to mesh together. These effects 
are due to the aging of the frost and are aptly included in the 
empirical relationship for $ . 

On the other hand, in Yamakawa, et al . (Reference 13) 
and Nakamura (Reference 14) data where T^ = 24 5° K and c{> a - 4 0% 
to 90%, the much higher wall temperature coupled with higher 
ambient humidities are expected to cause the formation of smaller 
water droplets which would not freeze instantaneously on the wall. 
Therefore, initially, the value of 3 for their data should be 
lower than that of Brian, et al. data (Reference 1), as the presence 
of the spherical ice particles in the frost layer is less probable. 

In this case, the water droplets tend to remain at the same size 
and due to the higher wall temperature, nucleation within the 
frost layer is not probable, even though the frost surface tempera- 
ture increases with time. Due to the foregoing reasoning, it is 
apparent that the expression for £ c should also include 
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the effect of wall temperature and the effect of the aging of 
frost. On the basis of the reasons given above, and for 
preliminary studies on frost formation, the expression for $ is 
proposed as 


6 C = 6 


•j . 9 + . 1 + 


2 2 
8CT - 


( . 526T + 315) 

w 


; 


(48) 


When the proposed Equation (48) is used for cal- 
culating the average thermal conductivity of the frost at T 
around 250°K, the results show a good agreement with the data 
by Yamakawa, et al . (Reference 13) and Nakamura (Reference 14). 
It may also be noted that at T = 80 °K, Equation (48) reduces 
to Equation (45) which gives excellent agreement for the data 
by Brian, et al . (Reference 1). Thus the expression for thermal 
conductivity with 3 C given by Equation (48) is used for the 
present frost formation model and the results are discussed in 
the last section of this report. 
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SECTION 3 


HEAT AND MASS TRANSFER COEFFICIENTS 

The diffusion of water vapor from the moist air across the 
frost surface and then into the interior of the frost layer is 
the mechanism by which the frost layer grows. The rate of growth 
will be determined by the rate of this mass transfer and the 
rate of heat transfer to and from the frost layer. To quantita- 
tively predict the growth of the frost, the rates of mass and 
heat transfer across the frost surface plane shown in Figure 1, 
must be determined. As in other heat and mass transfer 
problems it is convenient to represent these convective transfer 
rates as a heat or mass transfer coefficient multiplied by a 
suitable potential difference, either temperature, T, or concen- 
tration, co , as appropriate. That is, q^ = h„(T -T ) , and 

O n a s 

rfV = h (to -to ) . This section deals with the several methods of 
t mas 

representing the heat and mass transfer coefficients that have 
appeared in the literature. These coefficients are analyzed and 
compared with the experimental data available. The poor com- 
parisons obtained indicated that modified coefficients were needed 
and consequently are developed in this section. 

Because of the limited availability of data, the heat and 
mass transfer coefficients could only be compared at two extreme 
experimental conditions. The first condition is the frost for- 
mation on a vertical wall under laminar natural convection. The 
other condition is the frost formation in a horizontal duct under 
turbulent forced convection. Therefore, only the heat and mass 
transfer coefficients for these conditions are examined. Heat 
and mass transfer processes for other conditions, such as the 
airfoil geometry, are not examined because no experimental data 
is available for comparison. By analyzing for the two extreme 
experimental conditions two tasks can be accomplished. The first 
is that the heat and mass transfer processes for the other con- 
ditions can be inferred, and secondly, the frost formation model 
can be severely tested for accuracy under the extreme conditions . 
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Because accurate transfer coefficients for the frost structure 
are not available in the literature some investigators have tried 
to use empirical correlations such as are used in other convective 
heat transfer problems. Two of these for turbulent forced 
convection in a duct are the Colburn equation (Reference 16) , 

StPr 2//3 - 0.023 Re -0 * 2 

or Nu = 0.023 Re^ * 3 Pr^ 3 = h^L/k , (49) 

n a 

where Nu, Pr, Re and St stand respectively for the non-dimension 
Nusselt, Prandtl , Reynolds and Stanton number, and the Dittus - 
Boelter equation (Reference 16) 

Nu = 0.023 Re 0 * 8 Pr° " 4 . (50) 

Once the heat transfer coefficient is known, some methods are 
available to find the mass transfer coefficients under forced 
convection (Reference 6) , 

V h m = C p ' (51 > 

and the Chilton - Colburn analogy (Reference 6) 

h„/h = CLe" 2/3 , (52) 

H m p 

where C is specific heat of air (J/g °K) , 

P 2 

h„ is heat transfer coefficient (w/m °C) , 

2 

h^ is mass transfer coefficient (g/m s) , and 

Le is Lewis number = p C D/k 

a p a 

These relationships have been slightly modified for special 
situations. Kays and Perkins (Reference 16) show that for the 
Prandtl number between 0.5 and 1.0 and a constant wall temperature, 
the heat transfer coefficient for turbulent forced convection in a 
duct can be better correlated by the equation 

Nu = 0.021 Re° * 8 Pr° * 6 , (53) 
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than by the Colburn equation. For a vertical plate under laminar 

natural convection, the Nusselt number, Nu* is correlated to the 

' z 

local Grashof number, Gr z by the equation (Reference 6) 

Nu = 0.421 Gr 1//4 Pr 1/2 . (54) 

z z 

To show that these transfer coefficients need modification, 
the heat transfer coefficient experimentally obtained by Yamakawa, 
et al . (Reference 13) and Nakamura (Reference 14) were compared 
with Equations (53) and (54), respectively. Equation (53) when 
compared to Yamakawa data gave errors biased by about 200%, while 
Equation (54) compared to Nakamura data gave about a 25% error bias. 
Further significance of these errors will also be shown in the 
next section. Other experimentors (References 1, 2 and 4) have 
attempted to measure the accuracy of the heat and mass transfer 
processes for frost. Unfortunately, they have not succeeded for 
all the possible experimental ranges in correlating the experimental 
data by empirical relations. 

To explain the experimental disagreements noted earlier and 
the difficulties of several experimenters, three postulates are 
advanced. The first is that the suction or blowing of the water 
vapor at the frost surface will affect the thermal boundary layer 
and thus also the heat transfer coefficient. Secondly, fog for- 
mation or nucleation in the boundary layer will enhance the water 
vapor mass transfer and possibly also enhance the heat transfer. 
Lastly, due to frost porosity, the heat and mass transfer co- 
efficients are affected by the effective area of the frost surface 
which equivalently can be termed as frost roughness. With these 
three postulates in hand, the results of Nakamura (Reference 14) 
and Yamakawa, et al. (Reference 13) , are analyzed to determine 
which of the postulates are significant for a particular experi- 
mental condition. Thus the modification of the heat and mass 
transfer coefficients can be quantitatively described for the 
two extreme experimental conditions and qualitatively described 
for all other experimental conditions. With the modified co- 
efficients in hand, the heat and mass fluxes to the frost surface 
can be calculated. 
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3.1 LAMINAR NATURAL CONVECTION 

The results from Nakamura {Reference 14) show that in laminar 
natural convection on vertical plates the surface roughness has 
no effect on the heat transfer coefficient because the sensible 
heat flux from the air to frost is mostly conductive, rather 
than convective. Also, the nucleation in the boundary layer has 
only a small effect bn the heat and mass transfer coefficients, 
because the amount of vapor mass flux is usually so small that there 
is only minute fog formation. However, for natural convection, the 
momentum, heat and mass transfer mechanisms are coupled together 
in the boundary layer. The analysis for this problem is taken from 
Okino and Tajima (Reference 17) for a vertical plate under laminar 
natural convection. After correcting an algebraic error* in their 


analysis, the results are as follows. 

h H H 


where 


Nu ir 


n = 


Gr H = 


V 


Sc- $ 


= n/5 

Gr H 


sh H = 


h m H 


P D 
a 


= n(l+w g ) 


240 ( 


20 

21 


+ Sc$ ) 


1/4 


g«- 


5 | T a -T s 


CO -CO 

a s 


. 6453 + 2.6453co 


-/ 


$ 2 + 2 — 
v z Pr 

1 + co 


0 + 1 — $ 


$ = 


1 - CO 


(55) , (56) 


(57) 


(58) 


(59) 


(60 


v is kinematic viscosity, 
g is gravitational acceleration (m/s 2 ), 

H is height of plate (m) , 
oo s is saturated absolute humidity, and 
Sc is Schmidt number. 

The Nusselt and Sherwood numbers (Nu 
equations fitted Nakamura's data (Reference 14) quite well. 


TT and Sh„) by these 

n ri 


*The algebraic error is that the term /Sc in n was placed in the 
denominator instead of the numerator of equation (57) above. 
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3.2 TURBULENT FORCED CONVECTION 


In the case of forced convection however, the mass transfer 
has only a minor effect on the heat transfer coefficient. This 
is shown as follows. A conservative estimate of the blowing 
parameter based on Yamakawa, et al . data (Reference 13) gives 

F = m t = .00002 (61) 



which, according to a graph presented by Rohsenow and Hartnett 
(Reference 6) , is so small that the exponent of the Reynolds 
number in Equation (53) remains constant at 0.8. A cursory look 
at Yamakawa' s Nusselt number (Reference 13) versus Reynolds' 
number plot in Figure 9 verifies this conclusion, which is in 
opposition to a conclusion drawn by Yamakawa, that the mass trans- 
fer affects the heat transfer coefficient. Yamakawa, et al . 
(Reference 13) thought that the mass transfer could be eliminated 
by making the inlet humidity the same as that at the frost surface 
and to establish their conclusion, they took the convective heat 
as the only heat source at the frost surface. The heat flux, 
however, as measured at the wall below the frost layer, is 
greater than the convective heat because, due to the temperature 
gradient driving force, the latent heat is also being released in 
the frost layer. There is no evidence that any correction was 
made for the latent heat contribution and thus, an incorrect 
interpretation of the data was obtained. Yamakawa (Reference 13) , 
also makes the point that the water vapor nucleates in the 
thermal boundary layer above the frost surface. This process 
releases the latent heat and causes a change in the temperature 
profile in the thermal boundary layer, as observed by Yamakawa, 
which would change the way the heat transfer coefficient is 
related to the Reynolds' number. However, in the turbulent forced 

4 

convection limit, (i.e.. Re 10 ) the heat and mass transfer 
mechanisms are effectively decoupled because the sensible heat 
flux is much greater than the heat flux generated by the release 
of the latent heat through freezing. Therefore, it can be assumed 
that the latent heat release rate has a negligible effect on the 
heat transfer coefficient. 
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On the other hand the heat transfer coefficient is affected 
by the frost surface roughness. From our analysis of Yamakawa, 
et al. data (Reference 13), the ratio of the effective area of 
the frost surface to a smooth area is found to be a constant 
at 1.95. The frost roughness variety has no effect on h fi because 
the boundary layer is already turbulent and h R is experimentally 
independent of the frost porosity. Therefore, the usual heat 
transfer coefficient is increased by a factor, r = 1.95, as is 
shown by the correlation, 

Nu = h H D e = 0.021 Re‘ 8 Pr' 6 with h* = h H (1.95) , (62) 

k 

a 

where h-j = heat transfer coefficient of empirical equation and 

ri 

h* = experimental heat transfer coefficient. 

Equation (62) is believed to be the most accurate to date 
for calculating the heat transfer coefficient over frost and is 
being used as the correlation for turbulent forced convection. 
There is very little scatter when applying Equation (62) to 
Yamakawa' s data, (Reference 13) even for points corresponding to 
relative humidities from 56% to 97% as shown in Figure 9. 
Furthermore, the points corresponding to the different humidities 
show no preferred bias. This supports the assumption that the 
nucleation in the boundary layer has a negligible effect on the 
heat transfer coefficient in this experimental data. 

The mass transfer coefficient could be obtained from heat 
transfer coefficient using the Chilton - Colburn analogy (Ref- 
erence 6). However, Yamakawa's data (Reference 13) shows a 
strong scatter for the analogy of the heat and mass transfer. 

It is difficult to determine what causes the scatter of the h* 
vs h^ (experimental mass transfer coefficient) plot shown in 
Figure 10, because the points are not identified by the corres- 
ponding humidities and Reynolds number, yet certain trends can be 
perceived. For a mass transfer coefficient less than about 30 

g/m s, where h = h*, h* can be calculated as if from the 
3 m m m 

Chilton - Colburn analogy shown as the upper line in Figure 10, 
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Nusselt Number 



Reynolds Number, Re 

Figure 9. Relationship between Nusselt's number and Reynolds* 
number (Reference 13) 


*{Kcal/m -hr*AC°) 




( 63 ) 


bfi/ 1 " 95 = h w /h* = C Le~ 2 / 3 
H m p 


h* 

m 


This correlation seems to rule out the effect of the surface 
roughness on the experimental mass transfer coefficient, h^. 

But this correlation is tentative because in the turbulent flow 
the mass transfer coefficient should be affected by the surface 
roughness in the same way as the heat transfer coefficient is. 


2 

At h^ = 30 g/m s one can see a very strong scatter in the 
data, and it is biased in such a way that one can propose a mass 
flux enhancement due to nucleation in the thermal boundary layer. 
By using Epstein and Rosner's approach (Reference 11), it can be 
shown that for the wide range of wall temperatures and ambient 
humidities, the water vapor mass flux increases by a factor from 
1.0 to a maximum of 2.95 over the usual mass flux of condensation 
on the wall surface. That is. 


1 < h*/h_ < 2.95 . (64) 

— m m — 

2 

For the mass transfer coefficient greater than 45 g/m s, 
the Chilton - Colburn analogy between h* and h£ can be modified to 

h*/h* = 1,5 C p Le_2/3 / (65) 

shown as the lower line in Figure 10, which means that the 

enhancement factor is h*/h = 1.3. That is, if the Chilton - 

m m 

Colburn analogy is still valid for a rough surface and a 
turbulent flow. 


3.3 SUMMARY OF HEAT AND MASS TRANSFER COEFFICIENT EQUATIONS 

So far we have considered the extremes of the experimental 
data. At one extreme we have the laminar natural convection, and 
at the other extreme we have the turbulent forced convection. In 
the region of combined natural and forced convection not only the 
heat and mass transfer are coupled but also the nucleation in 
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the boundary layer and the frost surface area affect the heat 

transfer to a significant degree. The surface roughness will 

affect the heat transfer in two ways. The first is that the 

surface roughness will cause the transition point (located 

between laminar and turbulent zones) to occur closer to the 

leading edge. Secondly, in the laminar zone, the effective 

surface area for the heat transfer coefficient is unity relative 

to the smooth surface area, but the effective area is 1.95 times 

greater than the smooth surface area in the fully turbulent 

zone. Although the effects of nucleation in the thermal boundary 

layer upon the heat transfer was considered negligible in both 

extreme experimental data, it may become significant in the 

combined natural and forced convection when the latent heat 

release rate in the boundary layer may become the same order of 

magnitude as the sensible heat. Lastly, from our studies of the 

data, the perturbations in the mass transfer coefficients is 

caused by the fog formation or nucleation in the boundary layer 

and by the frost surface roughness. The maximum enhancement in the 

mass transfer coefficient appears to be at the ratio, h*/h =1.3. 

m m 

For the comparison of the frost formation model to experi- 
mental data, Equations (55) to (60) , (62) and (65) will be used 

tentatively. The frost formation model however, can be immediately 
generalized to other air flow geometry, and surface roughness 
conditions, if the heat and mass transfer coefficients can be 
obtained for these conditions. To obtain the transfer coeffi- 
cients for an arbitrary airfoil geometry, computer programs need 
to be developed to couple the boundary layer models with the 
potential flow models. There are such programs in the literature 
(References 18, and 19) but none has yet satisfied our requirements 
to include the effect of water vapor mass flux nucleation and 
surface roughness on the flow properties near the wall. Some 
literature (References 20, and 21) provides ways of calculating 
the influence of surface roughness and the water vapor mass flux 
on the transfer coefficients. Other literature (References 22 
and 23) provide ways of calculating the enhancement of the 
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mass transfer coefficient due to nucleation in the boundary layer. 
Only a small amount of literature is available for analyzing the 
heat and mass transfer problem for combined natural and forced 
convection. A synthesis of this literature will be required for 
the study of overnight frost formation on an aircraft. 
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SECTION 4 

THE SIMULATION OF FROST FORMATION 


4.1 FROST FORMATION MODEL 


In Section 2, the frost thermal conductivity was formulated 
so that the heat flux through the frost layer could be calculated. 

In Section 3, the heat and mass transfer coefficients were 
formulated so that the heat flux to the frost surface from the 
air could be calculated. By matching the two heat fluxes at the 
frost surface the temperature of the. frost surface can be calculated. 
The frost surface temperature will change with time because the 
frost layer is becoming denser and thicker, which in turn affects 
the computations of the heat flux in the frost layer. During the 
process of frost formation it is assumed that the part of the 
water vapor transported to the surface freezes at the surface so 
as to increase the thickness of the frost. The remaining water 
vapor is diffused into the existing frost layer before it freezes. 
Since the frost density is assumed to be spatially invariant in a 
direction normal to the plate, the water vapor diffusion flux 
entering the frost surface from the surrounding air is given by 


ih 


ds 


3p f 

dt X s 


( 66 ) 


where Pf and x g stand for water vapor diffusion mass flux at 

the surface, frost density and frost thickness respectively. 

Also, the water vapor mass flux, m^ will obey the diffusion 
Equation (18) and will be driven by the temperature gradient 
throughout the frost layer. 

The water diffusion flux is given by Equation (22) 



Equation (67) can be evaluated at the frost surface if the 

frost surface temperature T , the temperature qradient dT 

s dx 
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and the frost thickness x g are known. Combining Equations (66) and 
(67) allows the frost density to be solved for as a function of 
time. 

To obtain a value for frost surface temperature, T , the quasi- 
steady-state heat equation for the frost layer is to be solved. 

This is, the equation 

dT 

K -— ■ — = q_ (68) 

dx ° 

with K given by Equation (42) is to be solved with the boundary 
conditions 

T = T at x = x and T = T at x = 0 ' 
s s w 

where 

<Jo = h H (V T .> + " p vs )/p a 

+ za (T - T ) • 
a s 


This gives 


J 


s 

K dT = 


x a 
so 


(70) 


w 


The temperature gradient at the frost surface is given by 
dT 


dx 


K(T s ) 


(71) 


The thickness of the frost can be computed from the weight 


of the frost (p 


x. 


) which is directly related to the mass transfer 


coefficient and the water vapor density by 


d(p f x s ) 

dt 


= h ( p 
in K a 


= h m ( “a 


- p vs l/ p a 


“s ) ' 


(72) 


Since these equations cannot be solved analytically, numerical tech- 
niques are used to solve these equations to obtain the thickness 
and the density of the frost. 
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4.2 THE NUMERICAL SCHEME FOR FROST FORMATION MODEL 


The scheme begins by assuming a small initial value of frost 
density, p^q and frost thickness, x g Q. The initial frost density 
is set equal to that of the ambient air, i.e., p 
the initial thickness is calculated by 


fO 


= p ( 1 + w ) and 
a a 


sO 


= h* 
m 


(03 - 03 ) 

a w 


(A 


t o )//p fo 


(73) 


where is the absolute humidity at wall temperature and Atg is 
set at one second. The frost surface temperature, T is obtained 


The frost surface temperature, 
from the monotonic nonlinear equation, 

T_ 


/ 


w 


h ] 


R dT= x s | h H (T a -T s ) + h m L S <“ a -“s> + £ °<T a 4 - T s 4 )' 


]• 


(74) 


In order to solve this equation in an accurate and efficient manner, 
the frost thermal conductivity given by Equation (42) is first 
interpolated by a parabolic function in a temperature range AT and 
integrated to get the first term in Equation (74), (i.e.. At is 

often greater than T g - T w ) . Then T g in Equation (74) is solved 
for by an interval-halving-iterative technique so that the left term 
of Equation (74) is equal to the right terms. In this way successive 
evaluations of the complicated thermal conductivity expression (i.e. , 
Equation (42) ) and numerical integrations can be completely avoided 
when iterating to get T . The temperature gradient at the frost 


surface is obtained by 


dT 


dx 


r T 

f 3 

T 

W 


K dT 


(75) 


K (T ) 
s 


The value for the frost density is predicted for a later time t n+ ^ 
by using the equation. 


3P. 


_ DB 


at 


(l-X)T X 
s s 



- 1 


dT 


dx 


(76) 
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which is solved by using a second order Runge Kutta "initial" 
scheme, and by a second order predictor-corrector "takeover" 

scheme. Finally, the value for the frost thickness is predicted 

/ 

for the same time as the frost density, (i.e., at t + ^) 
by the equation. 


p f x s =\ n+1 h m (w a"“s )dt + p fO x sO ' 
*0 


( 77 ) 


where the rectangular rule is used for the integration, in 

acknowledgement of the fact that h and a* are slowly changing 

m s 

functions of time. At the 'same particular time, t a new 

frost surface temperature is calculated by returning back to 

Equation (74) incrementing n, and using the new values of frost 

density and thickness x . 

i s 
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SECTION 5 


COMPARISON OF THE MODEL WITH THE AVAILABLE EXPERIMENTAL DATA 

In this section, three separate classes of experimental data 
which are considered to be complete and the transfer coefficients 
considered to be well known are described and compared with re- 
sults predicted by the frost formation model. 

The first set of experimental data that the model predicts 
quite accurately is that of Brian, et al . (Reference 1) . Since 
the typical Reynolds number is greater than 10,000 and the boundary 
layer is turbulent in the duct in which their experiment was 
conducted, the correlation for the forced convection mentioned 
earlier by Equations (62) and (65) for the heat and mass transfer 
coefficients is used here. The roughness factor, r, is set at 
1.9 S'* and the enhancement factor is set at 1 . 3 . It may be noted 
that these coefficients are based on the experimental data of 
Yamakawa, et al. (Reference 13) . Table II lists the input data 
from Brian et al . (Reference 1) for the frost formation model and 
also the symbols identifying the experimental data in the next 
three figures. 

Results from the model have been compared with the experimen- 
tal data giving the frost weight, density, and thickness as 
functions of time under different experimental conditions. The 
results are given in Figures 11, 12 and 13 where the computer 
generated solid curves are plotted along with the experimental 
data. The figures show excellent agreement. 

When the standard correlations given by Equations (49) and 
(52) instead of the modified Equations (62) and (65) are used in 
this set of data, the corresponding results are shown by the 
dashed curves. The comparison of the solid and the dashed curves 
based on the same experimental data shows that poor correlations 
are obtained for standard heat and mass transfer coefficients. 

The second set of experimental data examined is that of 
Yamakawa et al. (Reference 13) for forced convection in a duct. 
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TABLE II 


Data Input to the Frost Formation Model for Comparison 
With Brian et al . Data for Forced Convection (Reference 1) 


Computer 
Curve No. 

Experimental 

Graph 

Symbol 

Reynolds 

Number 

Relative 

Humidity 

(%) 

Wall 

Temperature (°K) 

Air 

Temperature (°K) 

1 

O 

14500 

26.1 

79.83 

297 

2 

A 

9316 

26.2 

79.83 

297 

3 

□ 

5603 

26.0 

79.83 

297 

4 

• 

14750 

20.8 

79.83 

297 

5 

▲ 

9247 

15.8 

79.83 

297 

6 

■ 

5625 

16.2 

79.83 

297 



Weight (g/cm 



Figure 11. Weight versus Time for Brian et al. Data (Reference 1) 




Density (g/cc) 



Figure 12. Density versus 


time for Brian et al . Data (Reference 1) 




Thickness (cm) 



Figure 13. 


Thickness versus Time for Brian et al. Data (Reference 1) 


The roughness factor and the enhancement factor are again set 
at 1.95 and 1.3 respectively. 


Table III provides the input data into the frost formation 
model and Figures 14, 15 and 16 show the results of the computer 
generated solid curves along with the experimental data of frost 
weight, density, and thickness as functions of time. The dashed 
lines, generated using the standard correlations again, provide 
inferior results compared to the modified correlations. 

Curves 1 to 4 , in Figures 14 to 16, show a fair agreement 
with the experimental data. However; curve 5 shows a sharp 
disagreement with the data. The disagreement is probably due to 
the following complications. Since the relative humidity for 
this run is nearly 100%, the fog starts forming before the water 
vapor enters into the boundary layer. This means that the en- 
hancement factor becomes unity. Secondly, the water droplets 
reaching the surface are large, which in turn cause the frost 
structure to consist of a relatively large number of spherical 
ice-droplets. This alteration in the frost structure lowers the 
thermal conductivity to under the computed value. Third, the 
frost surface temperature quickly reaches the melting temperature, 
which causes the frost thickness to stabilize and the frost 
density to increase over their modeled values. The possibility 
of modifying the present model to accommodate these three factors 
will be considered for future work. 

The third set of experimental data examined is that of 
Nakamura, et al . (Reference 14) for natural convection on a 
vertical plate. The heat and mass transfer coefficients are given 
by Equations (55) and (56) , and the thermal conductivity is the 
same as that used for the Yamakawa data (Reference 13) . Table IV 
provides the input data into the frost formation model and 
Figures 17, 18 and 19 show the results of the computer generated 
curves along with the experimental data of frost weight, density, 
and thickness. Experimental data in the 5th and 6th entries of 
Table IV are not available for Figures 17 and 18. 
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TABLE III 


Data Input to the Frost Formation Model for Comparison With Yamakawa, 
et al. Data for Forced Convection in a Duct (Reference 13) 


Computer 
Curve No. 

Experimental 

Graph 

Symbol 

' 1 1 

Air 

Velocity 

(m/s) 

Relative 

Humidity 

(%) 

Wall j 

Temperature ( °K ) ; 

Air 

Temperature (°K) 

1 

□ 

1.6 

56.2 

251 

284 

2 

A 

5.0 

56.2 

251 

284 

3 

0 

5.0 

70.0 

251 

284 

4 

▲ 

7.3 

56.2 

251 

284 

5 

* 

7.3 

97.4 

251 

284 


Note: To calculate Reynolds Number, the hydraulic diameter is taken as 

D - 20(1. 7-x ) / ( 11 . 7-x ) cm. 
e s s 



Weight (g/cm 




Figure 14. Weight versus Time for Yamakawa, et al . Data (Reference. 13) 




Figure 15. Density 
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TABLE IV 


Data Input to the Frost Formation Model for Comparison With Nakamura 
Data for Natural Convection on Vertical Plate (Reference 14) 


Ul 


Computer 
Curve No. 

Experimental 

Graph 

Symbol 

Plate 

Height 

(meters) 

Ambient 
Relative 
Humidity (%) 

Wall 

Temperature (°K) 

Air 

Temperature (°K) 

1 

0 

0.085 

74 

262.46 

301.06 

2 

A 

0.085 

76 

259.66 

296.36 

3 

a 

0.085 

75 

257.36 

292.56 

4 

• 

0.085 

61 

260.76 

286.36 

5 

A 

0.040 

53 

256.56 

291.96 

6 

■ 

0.040 

38 

258.76 

291.66 
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Time (minutes) 

Figure 17. Weight versus Time for Nakamura Data (Reference 14) 
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Figure 18. Density versus Time for Nakamura Data (Reference 14) 




Thickness (cm) 



Figure 19. Thickness versus Time for Nakamura Data (Reference 14) 



It is interesting to note that the curves 2, 3 and 4 in Figure 17 
show a very good agreement with the experimental data. This 
verifies the accuracy of the mass transfer coefficient calculated 
in this experiment. However, the curves 1 and 2 in Figure 18 
show a disagreement with the experimental data of frost density 
and thickness versus time. The main reason seems to be that the 
frost surface starts melting very early, causing the thickness 
to stabilize sooner and the frost density to increase over its 
modeled value. 


SECTION 6 

SUMMARY AND RECOMMENDATIONS 

Keeping in view the complexity of the problem and the 
sparse experimental data available, the results obtained to date 
are very encouraging. The good agreement of the proposed model 
with the three sets of experimental data available has encouraged 
us to further develop this model for use in evaluating the aero- 
dynamic penalties on an aircraft having frost on its wings. 

Although the equations used in the present model for calculating 
the frost density as a function of time has done remarkably well 
for most frost formations, new equations must be developed to 
calculate the change in frost density due to liquid water seepage 
into the frost and for a non-spatially-uniform frost density. 

To realistically apply the frost formation model to an airfoil, 
additional considerations must also be given to the enthalpy 
rate terms for the aircraft skin and the frost, the complex 
black-body radiation to space, and the transfer coefficients for 
the various flow regimes and airfoil geometries. The formulation 
of the frost thermal conductivity used in the model has been 
placed on more theoretical grounds than in any previous 
approaches. As a consequence, the empiricism has been reduced to 
determining a structure parameter, $ , tentatively given by Equation 
48, which covers nearly all the experimental data examined to 
date. However, the validity of the structure parameter $ c 
needs to be verified at high frost density for the various 
values of wall-temperature. 

In summarizing, the present model for frost formation was 
constructed with a view towards its application to the overnight 
frost formation on an airfoil. The model contains many complexities 
not considered by any previous investigator. To obtain an 
improvement in the present model of frost formation would require 
more complete experimental data. In particular, specific 
experiments should be performed to improve our knowledge of the 
heat and mass transfer coefficients, improve the frost density 
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It is recommended that further development of the present 
frost model for utilization in the Aircraft Takeoff Simulation 
Program proceed along the following lines. Additional terms 
should be added to the heat balance equations to allow for 
radiation heat loss of the skin of the airfoil to space. This 
is the mechanism in nature that produces the rapid cooling of 
the airfoil, relative to the surrounding air, and results in 
frost formation. The resulting model could then be evaluated 
by comparison of computer results with experimental data. The 
experimental data, consisting of frost density and thickness 
values, would be obtained from measurements on a flat plate 
exposed to the environment during atmospheric conditions 
conducive to frost formation. Monitoring of the atmosphere 
during the frost formation period would provide the input data 
to test the computer model versus the actual observational data. 
The predictive capabilities of the model could also be tested. 
Using meteorological predications of the overnight temperature, 
pressure and cloud cover profiles , the computer model would be 
run to determine the anticipated frost accumulation by morning. 
This result would then be compared with the actual accumulation 
to assess the frost prediction capabilities of the model. An 
accurate frost prediction model giving reliable estimates of 
frost thickness and density would provide significant benefits 
to the aviation community. 

Finally after having established the validity of the frost 
model it should be incorporated into an aircraft simulation 
program to assess aerodynamic penalties during takeoff. In 
addition to providing the initial frost condition at takeoff, 
the model would be used throughout the takeoff phase of 
operation to determine the frost dissipation, or further 
accumulation, on the airfoil as the speed of the aircraft 
increases . 
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If 


rate modeling, especially with respect to water seepage, and 
obtain a more refined empirical formulation of the structure 

parameter, $ . 

c 
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